The FLUKA code: an overview 
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Abstract. FLUKA is a multipurpose MonteCarlo code which can transport a variety of 
particles over a wide energy range in complex geometries. The code is a joint project of INFN 
and CERN: part of its development is also supported by the University of Houston and NASA. 
FLUKA is successfully applied in several fields, including but not only, particle physics, cosmic 
ray physics, dosimetry, radioprotection, hadron therapy, space radiation, accelerator design and 
neutronics. The code is the standard tool used at CERN for dosimetry, radioprotection and 
beam-machine interaction studies. Here we give a glimpse into the code physics models with a 
particular emphasis to the hadronic and nuclear sector. 


1. Introduction 

Reliable calculation models are becoming an essential tool in nuclear physics applications. In 

particular there exist cases, concerning both fundamental and applied research, where accurate 

Monte Carlo simulations are required, fluka [1] is one of the existing codes to simulate 

transport and interaction of particles in matter which is capable of being a complete multipurpose 

tool, fluka is able to treat hadron-hadron, hadron-nucleus, nucleus-nucleus, neutrino, 
electromagnetic, and /i interactions up to 10000 TeV. Charged particle transport (handled in 
magnetic field too) includes all relevant processes. It also manages interaction and transport of 
neutrons down to thermal energies. 

In this review we shall give a description of the hadronic and nuclear physics models of fluka 

. Full descriptions of these models and of their applications can be found in the literature (see 
the web page, www.fluka.org). 



2. The nuclear models in FLUKA 

fluka is based, as far as possible, on original and well tested microscopic models. Due to this 
“microscopic” approach, each step is self-consistent and has solid physical bases. Performances 
are optimized by comparison with particle production data at single interaction level. No tuning 
whatsoever is performed on “integral” data, such as calorimeter resolutions, thick target yields, 
etc. Therefore, final predictions are obtained with a minimal set of free parameters, fixed for all 
energies and target/projectile combinations. Results in complex cases as well as scaling laws and 
properties come forth naturally from the underlying physical models and the basic conservation 
laws are fulfilled a priori. 

The basic building block is the description of the hadron-nucleon (h-N) interaction over a 
wide energy range. This is essential to achieve a sound description of the hadron-nucleus and 
nucleus-nucleus interaction. 

3. Hadron— nucleon interaction models 

Elastic, charge exchange and strangeness exchange reactions are described as far as possible 
by phase-shift analysis and/or fits of experimental differential data. Standard eikonal 
approximations are often used at high energies. 

At the low energy end (below 100 MeV) the p-p and p-n cross sections rapidly increase with 
decreasing energy. The n-p and the p-p cross sections differ by about a factor three at the 
lowest energies, as expected on the basis of symmetry and isospin considerations, while at high 
energies they tend to be equal. 

The total cross section for the two isospin components present in the nucleon-nucleon 
amplitude is given by: 
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The same decomposition can be shown to apply also for the elastic and the reaction cross 
sections. 

The cross section for pion-nucleon scattering is dominated by the existence of several direct 

resonances, the most prominent one being the A(1232). Given the pion isotopic spin ( T = 1), 

the three n charge states correspond to the three values of T z . Thus, in the pion-nucleon system 

two values of T are allowed : T = ^ and T = |, and two independent scattering amplitudes, 

Ai and A 3 , enter in the cross sections. Using Clebsch-Gordan coefficients all differential cross 
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sections can be derived from the three measured ones: a (ir + p — » tt + p), a(n~p — > vr“p), and 
a ( 7 T~p — > 7 r°n). 

As far as particle production is concerned (inelastic hadron-nucleon interactions), the 
description becomes immediately more complex. Two families of models are adopted, depending 
on the projectile energy: those based on individual resonance production and decays, which cover 
the energy range up to about 5 GeV, and those based on quark/parton string models, which 
can provide reliable results up to several tens of TeV. 

3.1. h-N interactions at intermediate energies 

The inelastic channel with the lowest threshold, single pion production, opens already around 
290 MeV in nucleon-nucleon interactions, and becomes important above 700 MeV. In pion- 
nucleon interactions the production threshold is as low as 170 MeV. Both reactions are normally 
described in the framework of the isobar model: all reactions proceed through an intermediate 
state containing at least one resonance. There are two main classes of reactions, those which 
form a resonant intermediate state (possible in 7r-nucleon reactions) and those which contain two 



particles in the intermediate state. The former exhibit bumps in the cross sections corresponding 
to the energy of the formed resonance. Examples are reported below: 

Ali + IV 2 -> N[ + A(1232) ^N^ + N^ + tt 
n + N -»• A(1600) -> it' + A(1232) -> 

— + tt' + it" + N' 

m + N 2 -> Ai(1232) +A 2 (1232) -> 

— > Ni + 7Tl + N 2 T 7T2 

Partial cross sections can be obtained from one-boson exchange theories and/or folding of 
Breit-Wigner with matrix elements fixed by N-N scattering or experimental data. Resonance 
energies, widths, cross sections, and branching ratios are extracted from data and conservation 
laws, whenever possible, making explicit use of spin and isospin relations. They can also be 
inferred from inclusive cross sections when needed. 

3.2. Inelastic h-N interactions at high energies 

As soon as the projectile energy exceeds a few GeV, the description in terms of resonance 
production and decay becomes more and more difficult. The number of resonances which should 
be considered grows exponentially and their properties are often poorly known. Furthermore, 
the assumption of one or two resonance creation is unable to reproduce the experimental finding 
that most of the particle production at high energies occurs neither in the projectile nor in the 
target fragmentation region, but rather in the central region, for small values of Feynman x 
variable. Different models, based directly on quark degrees of freedom, must be introduced. 

The features of “soft” interactions (low-pr interactions) cannot be derived from the QCD 
Lagrangian, because the large value taken by the running coupling constant prevents the use 
of perturbation theory. Models based on interacting strings have emerged as a powerful tool 
in understanding QCD at the soft hadronic scale, that is in the non-perturbative regime. An 
interacting string theory naturally leads to a topological expansion. The Dual Parton Model 
(DPM) [2] is one of these models and it is built introducing partonic concepts into a topological 
expansion which explicitly incorporates the constraints of duality and unitarity, typical of 
Regge’s theory. In DPM hadrons are considered as open strings with quarks, antiquarks or 
diquarks sitting at the ends; mesons (colourless combination of a quark and an antiquark qq) 
are strings with their valence quark and antiquark at the ends. At sufficiently high energies 
the leading term in the interactions corresponds to a Pomeron (IP) exchange (a colourless 
closed string exchange), which has a cylinder topology. When an unitarity cut is applied to 
the cylindrical Pomeron, two hadronic chains are left as the sources of particle production, as 
pictorially represented in Fig.l. 

As a consequence of colour exchange in the interaction, each colliding hadron splits into two 
coloured system, one carrying colour charge c and the other c. The system with colour charge 
c (c) of one hadron combines with the system of complementary colour of the other hadron, 
to form two colour neutral chains. These chains appear as two back-to-back jets in their own 
centre-of-mass systems. 

The exact way of building up these chains depends on the nature of the projectile-target 
combination (baryon-baryon, meson-baryon, antibaryon-baryon, nreson-nreson). Further 
details can be found in the original DPM references [2] or in [3]. 

The chains produced in an interaction are then hadronized. DPM gives no prescriptions 
on this stage of the reaction. All the available chain hadronization models, however, rely on 
the same basic assumptions, the most important one being chain universality. This implies 
that chain hadronization does not depend on the particular process which originated the chain, 
and until the chain energy is much larger than the mass of the hadrons to be produced, the 
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Figure 1. Left: cylindrical topology of pomeron exchange. Right: the two chains resulting 
from the unitarity cut for a p-p collision. 


fragmentation functions (which describe the momentum fraction carried by each hadron) are 
the same. As a consequence, fragmentation functions can in principle be derived from hard 
processes and e + e _ data and the same functions and (few) parameters should be valid for all 
reactions and energies. Actually mass and threshold effects are non-negligible at the typical chain 
energies involved in hadron-nucleus reactions. As a last step in Dual Parton Model, a transverse 
momentum has to be added and this is done starting from uncertainty considerations. An 
example of fluka performances is given in Fig. 2, where the double differential cross section for 
inclusive A production in K-p collision at 10 GeV incident laboratory energy is reported. 

4. Main steps of hadron— nucleus interactions 

High energy hadron-nucleus (h-A) interactions can be schematically described as a sequence of 
the following steps: 

• Glauber-Gribov cascade 

• (Generalized) IntraNuclear cascade ((G)INC) 

• Preequilibrium emission 

• Evaporation/Fragmentation/Fission and final deexcitation 
4-1- The Glauber-Gribov cascade and the formation zone 

The Glauber formalism [4, 5] provides a powerful and elegant method to derive elastic, quasi- 
elastic and absorption h-A cross sections from the free h-N cross section and the nuclear ground 
state only. Inelastic interactions are equivalent to multiple interactions of the projectile with u 
target nucleons. The number of such “primary” interactions follows a binomial distribution (at 
a given impact parameter, b ): 


Pruib) = [^P r v (b)[\-P r (b)} A - V 

where P r (b) = ahNrT r (b), and T r (b) is the profile function (folding of nuclear density and 
scattering profiles along the path). On average: 
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Figure 2. Comparison of simulated vs. experimental double differential cross section for 
inclusive A production in K-p collision at 10 GeV incident laboratory energy. 
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The Glauber-Gribov model [6, 7, 8] represents the diagram interpretation of the Glauber 
cascade. The v interactions of the projectile originate 2u chains, out of which 2 chains (valence- 
valence chains) struck between the projectile and target valence (di)quarks, 2{y — 1) chains 
(sea-valence chains) between projectile sea q — q and target valence (di)quarks. 

The distribution of the projectile energy among many chains naturally softens the energy 
distributions of reaction products and boosts the multiplicity with respect to hadron-hadron 
interactions. In this way, the model accounts for the major A-dependent features without any 
degree of freedom, except in the treatment of mass effects at low energies. 

The Fermi motion of the target nucleons must be included to obtain the correct kinematics, in 
particular the smearing of p j- distributions. All nuclear effects on the secondaries are accounted 
for by the subsequent (G)INC. 

The formation zone concept is essential to understand the observed reduction of the re- 
interaction probability with respect of the naive free cross section assumption. It can be 
understood as a “materialization” time. At high energies, the “fast” particles (in the emulsion 
language particles with /? >0.7) produced in the Glauber cascade have a high probability to 
materialize already outside the nucleus without triggering a secondary cascade. Only a small 
fraction of the projectile energy is thus left available for the INC and the evaporation. 

The Glauber cascade and the formation zone act together in reaching a regime where the 
“slow” part of the interaction is almost independent of the particle energy. The average 
multiplicity and its variance are directly related to the distribution of primary collisions as 


predicted by the Glauber approach. Due to the very slow variation of h-N cross section from 
a few GeV up to a few TeV, the Glauber cascade is almost energy independent and the rise in 
the multiplicity of “fast” particles is related only to the increased multiplicity of the elementary 
h-N interactions. 

Due to the onset of formation zone effects, most of the hadrons produced in the primary 
collisions escape from the nucleus without further reinteractions. Further cascading only involves 
the slow fragments produced in the target fragmentation region of each primary interaction, and 
therefore it tends quickly to saturate with energy as the Glauber cascade reaches its asymptotic 
regime. This trend is well reflected in the average multiplicity (and multiplicity distribution) 
of “gray” tracks (charged particles with 0.3 < /3 < 0.7), which are mostly protons produced in 
secondary collisions during the INC and preequilibrium phases. 

At the end of the cascading process, the residual excitation energy is directly related to the 
number of primary and secondary collisions which have taken place. Each collision is indeed 
leaving a “hole” in the Fermi sea which carries an excitation energy related to its depth in 
the Fermi sea. Evaporation products, as well as residual excitation functions, should reach an 
almost constant condition as soon as the Glauber mechanism and the formation zone are fully 
developed. This can actually be verified by looking at the production of “black” tracks (charged 
particles with (3 < 0.3), which are mostly evaporation products. 

4-2. (Generalized) IntraNuclear cascade 

At energies high enough to consider coherent effects as corrections, a h-A reaction can be 
described as a cascade of two-body interactions, concerning the projectile and the reaction 
products. This is the mechanism called IntraNuclear Cascade (INC). INC models were already 
developed at the infancy of the computer era with great success in describing the basic features 
of nuclear interactions in the 0.2-2 GeV range. Modern INC models had to incorporate many 
more ideas and effects in order to describe in a reasonable way reactions at higher and lower 
energies. Despite particle trajectories are described classically, many quantistic effects have to 
be incorporated in these (G)INC models, like Pauli blocking, formation time, coherence length, 
nucleon antisymnretrization, hard core nucleon correlations. A thorough description of the 
(G)INC model used in fluka can be found in [1, 3]. 

4-3. Preequilibrium 

At energies lower than the it production threshold a variety of preequilibrium models have been 
developed [9] following two leading approaches: the quantum-mechanical multistep model and 
the exciton model. The former has a very good theoretical background but is quite complex, 
while the latter relies on statistical assumptions, and it is simple and fast. Exciton-based models 
are often used in Monte Carlo codes to link the INC stage of the reaction to the equilibrium 
one. 

In the fluka implementation the INC step goes on until all nucleons are below a smooth 
threshold around 50 MeV, and all particles but nucleons (typically pions) have been emitted or 
absorbed. At the end of the INC stage a few particles may have been emitted and the input 
configuration for the preequilibrium stage is characterized by the total number of protons and 
neutrons, by the number of particle-like excitons (nucleons excited above the Fermi level), and 
of hole-like excitons (holes created in the Fermi sea by the INC interactions), and by the nuclear 
excitation energy and momentum. All the above quantities can be derived by properly recording 
what occurred during the INC stage. The exciton formalism of fluka follows that of M. Blann 
and coworkers [10, 11, 12, 13], with some modifications detailed in [3]. 



4-4 ■ Evaporation, fission and nuclear break-up 

At the end of the reaction chain, the nucleus is a thermally equilibrated system, characterized 
by its excitation energy. This system can “evaporate” nucleons, fragments, or 7 rays, or can 
even fission, to dissipate the residual excitation. 

Neutron emission is favoured over charged particle emission, due to the Coulomb barrier, 
expecially for medium-heavy nuclei. Moreover, the excitation energy is higher in heavier nuclei 
due to the larger cascading chances and to the larger number of primary collisions in the Glauber 
cascade at high energies. The level density parameter a ~ A / 8 MeV is higher too, thus the 
average neutron energy is smaller. Therefore, the neutron multiplicity is higher for heavy nuclei 
than for light ones. 

The fluka evaporation module is based on the standard Weisskopf-Ewing formalism [14]. 
Latest improvements [1] are represented by: i) adopted state density expression 

p oc exp (2 \ZaU)/U^ (where U is the emitting nucleus excitation energy), ii) no Maxwellian 
approximation for energy sampling, iii) competition with 7 emission, iv) sub-barrier emission. 
Neutron and proton production are marginally affected, while residual nuclei production and 
alpha emission are nicely improved. 

For light residual nuclei, where the excitation energy may overwhelm the total binding energy, 
a statistical fragmentation (Fermi Break-up) model is more appropriate (see [1, 3, 15] for the 
fluka implementation) . 

The evaporation/fission/break-up processes represent the last stage of a nuclear interaction 
and are responsible for the exact nature of the residuals left after the interactions. However, 
for a coherent self-consistent model, the mass spectrum of residuals is highly constrained by the 
excitation energy distribution found in the slow stages, which in turn is directly related to the 
amount of primary collisions and following cascading which has taken place in the fast stages. 

4-5. Ion- Ion interactions in fluka 

The fluka implementation of suitable models for heavy ion nuclear interactions has reached 
an operational stage. At medium/high energy (above a few GeV/n) the dpmjet model is 
used, dpmjet [16] is a Monte Carlo model for sampling hadron-hadron, hadron-nucleus and 
nucleus-nucleus collisions at accelerator and cosmic ray energies (E^ from 5-10 GeV/n up 
to 10 11 GeV/n). dpmjet is built upon the same principles used for the high energy part of 
fluka , i.e. is based on the two component Dual Parton Model in connection with the Glauber 
formalism, fluka implements both dpmjet-II.53 and dpmjet-III[17] as event generators to 
simulate nucleus-nucleus interactions. A comparison of data vs simulation is shown in Fig. 3. 
De-excitation and evaporation of the excited residual nuclei is performed by calling the FLUKA 
evaporation module. 

The DPMJETnrodel is not valid for energies below a few GeV/nucleon. For this reason, an 
interface to the rqmd-2.4 model was developed to enable fluka to treat ion interactions from 
rTOO MeV/n up to 5 GeV/n The rqmd-2.4 [18] is a relativistic model based on “Quantum 
Molecular Dynamics” (QMD). This is an approach where individual nucleons evolve according 
to an effective Hamiltonian, involving two- and three-body interaction terms, rqmd-2.4 has 
been successfully applied by the original authors to relativistic AA particle production over a 
wide energy range, from ~ 0.1 GeV/n up to several hundreds of GeV/n. Several important 
modifications have been implemented in the RQMD code, in order to ensure energy-momentum 
conservation taking into account experimental binding energies, and to provide meaningful 
excitation energies for the residual fragments. The results of this modified model in the energy 
range of interest for the applications described here can be found in [19] . 
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Figure 3. Comparison of simulated vs. experimental pseudorapidity distributions of charged 
hadrons produced in Au-Au collisions at a c.nr. energy of 130 GeV/A for different values of the 
“centrality” parameter of the collisions. Data points (labels marked with B) are taken from the 
Brahms, Phobos and Phenix experiments at RHIC. 


5. Developments in progress 

A new QMD code is also being developed from scratch to describe A-A collisions taking into 
account both the effect of stochastic scattering between all the involved nucleons and the effect of 
the nuclear potential which acts on each nucleon. Protons and neutrons are described as gaussian 
wave packets of fixed widths; the total nuclear wave function is approximated by the product of 
single nucleon wave functions. Regarding the Hamiltonian, each group of QMD model developers 
makes its particular choices (see, for instance, refs. [20, 21, 22, 23, 24, 25]); our starting point 
is a non-relativistic phenomenological potential, based on Skyrme interaction, supplemented 
by surface and symmetry terms; we also add the electromagnetic repulsion between protons, 
crucial to determine low-energy nuclear trajectories. In principle one can build a Hamiltonian 
as sophisticated as desired to improve the nuclear physical description; in practice however one 
has to meet CPU time requirements. 

The parameters of the model are fixed in order to reproduce as accurately as possible the 
observed nuclear ground state properties. We emphasize that this result can be achieved only 
approxinratively, because one has to describe with only a few parameters a great variety of 
nuclei, ranging from the lightest to the heaviest ones. 

We underline that QMD codes are based on Monte Carlo simulations: the cross-sections for 
A-A interactions are obtained as mean values from hundreds and hundreds of events, each of 


which should involve different starting nuclear configurations. In practice, it turns out that one 
can simulate a wide variety of different events with only a few initial configurations, randomly 
rotated. The choice and the storage of reasonable initial configurations can be pre-computed. 

We are now completing the coupling of the dynamical nuclear evolution predicted by 
our QMD (which gives the description of the first stage of the reaction) with the FLUKA 
preequilibrium module, to describe the deexcitation of the fragments formed and to study deeply 
the fragmentation process, implementing suitable models. Preliminary results are reported in 

[27]- 

Moreover, a promising task is represented by the coupling of FLUKA with a Monte Carlo 
code [28] developed at Milan University and based on Boltzmann Master Equation theory, as a 
tool to treat ion-ion interactions below 100 MeV/n. An example of the performance is reported 
in Fig. 4. Other details about this development can be found in [30]. 



Figure 4. Comparison of the prediction of BME Model of ref. [28] for double differential neutron 
yield resulting from 20 Ne collisions against 165 Ho at a total incident laboratory energy of 292 
MeV. Data are from ref. [29]. 
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